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^ ■ ABSTRACT 

"*^ I We investigate numerically the hydrodynamic instability of an ionization front 

^ I (IF) accelerating into a molecular cloud, with imposed initial perturbations of 

different amplitudes. When the initial amplitude is small, the imposed perturba- 
. _ tion is completely stabilized and does not grow. When the initial perturbation 

IT) ' amplitude is large enough, roughly the ratio of the initial amplitude to wave- 

in 



length is greater than 0.02, portions of the IF temporarily separate from the 

■^ I molecular cloud surface, locally decreasing the ablation pressure. This causes 

\^ [ the appearance of a large, warm HI region and triggers nonlinear dynamics of 

O I the IF. The local difference of the ablation pressure and acceleration enhances 

'^ I the appearance and growth of a multimode perturbation. The stabilization usu- 

I I ally seen at the IF in the linear regimes does not work due to the mismatch of 

H I the modes of the perturbations at the cloud surface and in density in HII re- 

52 I gion above the cloud surface. Molecular pillars are observed in the late stages of 

the large amplitude perturbation case. The velocity gradient in the pillars is in 

KJ reasonably good agreement with that observed in the Eagle Nebula. The initial 

Vj ' perturbation is imposed in three different ways: in density, in incident photon 

number flux, and in the surface shape. All cases show both stabilization for a 

small initial perturbation and large growth of the second harmonic by increasing 

amplitude of the initial perturbation above a critical value. 
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1. INTRODUCTION 

The shapes of the surfaces between molecular clouds and the HII regions around massive 
stars typically feature elongated structures, commonly referred to as pillars, columns, spikes, 
or elephant trunks (Hester et al. 1996; Pound 1998; Pound et al. 2003; Bally 2003). The 
surface is a photoionization front driven by the strong UV radiation from the OB stars. 
One well-known example is the Eagle Nebula which has three large, molecular pillars near a 
small group of O stars (Hillenbrand, Massey, Strom, & Merrill 1993). Although a number of 
theoretical and numerical studies have been done, the formation mechanism is still not fully 
understood. 

The OB stars are hot, with photospheres at temperature of a few tens of thousands of 
Kelvins (Dors & Copetti 2003), and give off high intensity UV photons. These UV photons 
irradiate the molecular cloud which surrounds the OB stars and photoevaporation occurs, 
resulting in a stratified structure. The photo evaporative fiow (ablated plasma) velocity is 
normal to the molecular surface (Hester et al. 1996). The region between the OB stars 
and the molecular cloud surface is the HII region in which the hydrogen gas is almost fully 
ionized and the photoionization and recombination to neutral atomic hydrogen occurs in 
steady state. The IF is a very thin layer because of the short mean free path for the incident 
photons on the cloud surface. At this cloud surface all the photons above the Lyman limit 
are absorbed. The photons below the ionization limit but above 11.2 eV cannot ionize the 
hydrogen atom in the ground state but can penetrate the IF and dissociate the molecular 
hydrogen in the underlying layers. The thickness of the dissociation front differs with each 
cloud. The last layer in the stratified structure is the molecular gas. Because of strong 
radiative cooling, the molecular cloud temperature is typically a few tens of Kelvin (Neufeld, 
Lepp, & Melnick 1995). The dynamics of such irradiated molecular clouds is thought to 
play a role in star formation (Bally 2003; Sugitani et al. 2002; McCaughrean & Andersen 
2002). Consequently, the outfiow dynamics from the molecular cloud into the HII regions is 
of considerable, general interest (Pound et al. 2003; Mizuta, Yamada, & Takabe 2002). 

The pillars in M16, popularly known as the Eagle Nebula, are at a distance of ~ 2 kpc 
((Hillenbrand, Massey, Strom, & Merrill 1993)). It is estimated that half of the ionizing 
radiation comes from a single 03-4 star and the rest mostly comes from the three other 
nearby stars (05-6) (Hester et al. 1996). The orientation of this system can be seen in Fig. 
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1 (taken from Pound (1998)). The total ionizing flux is estimated to be 5 ~ 1.2 x lO^*' s^^. 
The temperature in the HIT region near the molecular cloud surface is 9500 K and the ionized 
hydrogen number density is 5 x 10^ cm~^ (Levenson et al. 2000). The thickness of the 
photodissociation region (PDR) is relatively thin because the large hydrogen number density 
(n(H) ~ 10^ cm~^ ) provides high optical depth. As measured by CO line observations, the 
velocity gradients along the long axes of the pillars (from "head" to "tail") are between 
-20.7 and +6.7 km s^^pc^^, with an average magnitude of 8.3 km s^^pc^^. Furthermore, 
the pillars are not in the same plane in the sky; the differing signs of the velocity gradients 
indicate inclination towards (positive gradient) or away (negative gradient) from the observer. 
Predictions of the classical Rayleigh- Taylor (RT) instability theory with regard to the velocity 
gradient (Frieman 1954) are incompatible with the observed gradient, at least for constant 
acceleration (Pound 1998). Another formation mechanism should be considered. 

The formation mechanism of pillars beside massive stars has been discussed by a num- 
ber of authors. There are two classes of models for their formation. One is that the pillars 
are formed due to dense, pre-existing cores in the molecular cloud (e.g. Reipurth (1983); 
Bertoldi & McKee (1990); Lefioch & Lazareff (1994); Wilhams, Ward-Thompson, & Whit- 
worth (2001)). In a interesting variation on this idea, Williams, Ward-Thompson, & Whit- 
worth (2001) also simulated the result of a radiatively driven, pre-existing, short pillar with 
the same density as the background molecular cloud (their Case II). Even in this case, the 
pillar grows. 

Another possibility is that the pillars are caused by hydrodynamic instabilities. Spitzer 
(1954) proposed a model based on the RT instability occurring at the contact discontinuity 
between the dense molecular cloud and the lower density, hot, photoevaporated plasma. He 
suggested that the pillars in the Eagle Nebula were the result from the nonlinear stage of the 
RT instability. Frieman (1954) estimated the time scale of the instability to be less than 10® 
years. Pottasch (1958) showed observations of several nebulae and estimated the age of the 
pillars by the RT model, assuming small initial amplitude of the perturbations. Vandervoort 
(1962) found another type of instability at the IF in a non-accelerating frame, the so-called 
IF instability, which, in his analysis, was present in the case of non-normal incident radiation. 
He derived a dispersion relation for perturbations growing by this process. Axford (1964) 
extended it to include recombination, which plays a crucial role in the HII region. The 
recombination in the ionized gas works to smooth the surface when the wavelength of the 
perturbation is much larger than the recombination length. 

Williams, Ward-Thompson, & Whitworth (2001) showed robust development of pho- 
toionized pillars with an isothermal model, although they considered a semi-infinite cloud 
and a non-accelerating IF. Williams (2002) derived a dispersion relation for the IF with a 
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non-normal incident radiation field including the effect of recombination in the ionized gas. 
Ryutov et al. (2003) found unstable modes with non-normal radiation for an accelerating 
IF, but did not include recombination. 

Most studies have assumed a semi-infinite molecular cloud and non-accelerating IF. 
Instabilities at an accelerating IF were studied numerically by Mizuta et al. (2005, hereafter 
Paper I). They show that, for imposed perturbations with small initial amplitude, i.e. in 
the linear regime, there was no significant growth. Large growth of a classical RT instability 
is observed when recombination in the HII region is not included. We study here in more 
depth the dynamics of an accelerating IF. The existence of the acceleration at the IF is quite 
different than most previous models. This paper is organized as follows. The numerical 
method and conditions are described in Sec. 2. The results and discussions are given is Sec. 
3, and the conclusion is given in Sec 4. 



2. MODEL : ACCELERATING IF 

The same physics and computational method used in Paper I are included in this study. 
The energy balance and magnetic pressure to prevent radiative collapse for the molecular 
cloud are considered. The equations we numerically solve are: 
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n— + nu-Vf = an{l-f)\J\-aBn^f, (5) 

V-J = -an(l-/)|J|, (6) 

where p is mass density, p is pressure, / is unit tensor, u is the velocity vector, and e 
is the specific internal energy. The equations (l)-(3) are mass, momentum and energy 
conservation with energy sources. The energy source terms due to recombination in the 
ionized region, absorption of the UV radiation from OB stars, and cooling in the molecular 
gas are qreAuvAmoU respectively. These energy sources are evaluated as 

Qre = {nff pB k T, 
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Quv = Wan{l- f)\J\, 

Qmoi = n^noi X 10 erg cm s , 

T = K/A;)x[4e/(7/ + 5)], 

Umoi = n (1 - /)/2, 

pe = 3{nf)kT + 2.5nmoikT, 

where T is the temperature in Kelvins rrip is the proton mass, k is the Bohzmann constant, 
and n and rimoi are atomic hydrogen and molecular hydrogen number density. We ignore 
metal line cooling for simplicity, since its cooling power has the same dependence as that of 
recombination (proportional to ionized hydrogen number density). The dissociation heating 
in the neutral region and other radiative processes are also ignored for simplicity. The heating 
function is W = 1.73 x 10~^^ erg which corresponds to the average energy deposited into the 
gas per absorbed ionizing photon, and leads to produce an isothermal temperature of T = 10'^ 
K in the ionized gas as an equilibrium state of photoionization heating and recombination 
cooling. The case B recombination coefficients, which are summation of all recombination 
coefficients of hydrogen except the recombination to the ground state, are assumed constant 
as = 2.6 X lO-i^cm^s-i and (3b = l.25aB (at T = 10^ K from Hummer & Seaton (1963)), 
where f^B includes the thermal velocity dependence of the rates of recombination and free- 
free collisional cooling. The equation of state (Eq. 4) includes a magnetic pressure term for 
dense gas to prevent radiative collapse (Ryutov & Remington 2002), where / = rii/n is the 
ionization fraction, n^ is the ionized hydrogen volume density, p^^ and p^^ are constant values. 
The index 7m (= 4/3) is also constant. Numerically, we do not take into account the atomic 
hydrogen state, assuming the dissociation front is thin as in the Eagle nebula. Equations 
(5) and (6) describe the evolution of / and the transport of the incident radiation, where 
a = 6 X lO^^^cm^ is the photoionization cross-section of hydrogen, and J is the number flux 
of ionizing photons, i.e., photons cm^^s^^. Since we use the "on-the-spot approximation" 
(i.e., photons emitted in recombination to ground state are immediately reabsorbed, whereas 
the photons emitted in recombination to second or higher levels are assumed to escape from 
the system), we do not consider the diffusive photon emission and transport. 

The 2D computational domain is {x xy) = (0.46 pcx3 pc). Uniform grid points (Ax = 
Ay = 2.5 X 10~^ pc) are used. Periodic boundary conditions are employed (x = and 
X = 0.46 pc). The other boundary conditions are open boundary which means zero gradient. 
An incident photon flux (Jq = ~5 x lO^^y cm^^s^^) is imposed from one boundary {y = 3 
pc), where y is the unit vector of the y axis. The rays are parallel to the y axis. A cloud of 
initial thickness of a quarter pc, and density of n{H2) = 10^ cm~^ is set 0.5 pc away from 
the boundary where the incident photons come in. A very dilute gas n{H) = 10 cm~^ is 
imposed in the other regions. Initially, the gas at y > 2.25 pc is in pressure equilibrium. 
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and the gas at y < 2.5 pc is isothermal (40 K). The constant parameters in the equation of 
state, such as, p^^ and Pj^^, are the mass density and thermal pressure, of the initial molecular 
cloud, respectively. 

We impose the perturbation in the cloud in three different ways. The first is a density 
perturbation in a layer of 0.125 pc thickness along the surface of the cloud, with the form 
n = no(l — Acos[27ra;/0.46 pc]), where A is the amplitude of the perturbation and tiq = 
10^ cm""^. We assume a wavelength of A = 0.46 pc, and consider amplitudes A = 0.2 
(model D2), 0.3 (D3), 0.4 (D4) and 0.5 (D5). The second method is a 30% amplitude 
perturbation in the incident photon number flux of form J = Jo(l ~ 0.3cos[27ra;/0.46 pc]), 
starting at 98, 102, 106, or 108 kyr and ending at 110 kyr. These are models P098, P102, 
P106, P108, respectively. This is similar to the method used in Paper I and in Mizuta et al. 
(2005b). The third method is to impose surface perturbations on the initial cloud, according 
toi/ = 2.5-Ccos[27ra;/0.46pc]. The amplitude (C) is 3.8 x 10^^ pc (model S038), 9.0 x 10"^ 
pc (S090), 1.4 X 10-2 pc (S140), and 1.9 x 10"^ pc (S190). 



3. RESULTS and DISCUSSION 

3.1. ID Dynamics 

We briefly show the results without any perturbations, before showing the results with 
perturbations. The dynamics without any perturbation is very simple. When the incident 
photon flux is turned on, a shock propagates through the molecular cloud (compression 
phase). After this shock breaks out of the back side of the cloud, a rarefaction passes back 
through the shocked cloud. Then an acceleration phase begins at about t = 100 kyr, as 
the cloud moves as a unit. Figure 2 shows one dimensional profiles of the hydrogen number 
density at the early phase of the dynamics at t = 0, 50, 100, 110, and 150 kyr. In the figure, 
all three phases are shown, although the phase at around t = 100 kyr to t = 110 kyr in 
which the rarefaction passes back through the compressed cloud is too short to show. 

The hydrogen number density profiles of an initially semi-infinite cloud are shown in 
Fig. 3, both with and without magnetic pressure. The initial total pressure without the 
magnetic component is lower than with it. Without magnetic pressure, the number density 
of the compressed cloud becomes about 1.5 times higher, since there is less pressure support. 
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3.2. Small Initial Amplitude Cases 

Figure 4 shows the evolution of the amplitude, i.e., half the peak-to- valley of the per- 
turbation of the / = 0.5 contour in the y direction for each case, for all three methods of 
imposing the perturbation. It should be noted that this does not always correspond to the 
amplitude of the perturbation on the cloud surface because of the effect of "separation of 
the IF" discussed below. In cases with small initial perturbations (models D2, D3, P106, 
P108, and S038), the perturbation does not grow but rather oscillates as observed in Paper 
1. This stabilization is caused by the density profile in the Hll region. The ablated (pho- 
toevaporated) gas in the concave region of the incipient bubble concentrates (gets weakly 
"focused" ) in the region above the bubble vertex. This concentration of gas subsequently 
absorbs more of the incoming ionizing photons, reducing the photon flux reaching the IF 
region at the cloud surface near the bubble vertex, compared with the spike region. This 
is because the recombination rate depends on the square of the ionized hydrogen number 
density in the HII region. The net result is that the ablation pressure locally decreases in 
the bubble region compared to the spike region, which acts to smooth out the perturbed 
surface, if the perturbation is small. This is also shown in Fig. 5a with a 2D plot of number 
density late in time (t = 460 kyr) for a small initial density perturbation (model D2). The 
IF is quite smooth, namely, no growth of the perturbation has occurred. The other cases 
(D3, P106, P108, and S038) in which the perturbation does not grow behave similarly. 



3.3. Large Initial Amplitude Cases 

When the perturbation initial amplitude is larger, roughly when the ratio of the initial 
amplitude to wavelength is greater than 0.02, the perturbations grow, after the acceleration 
phase begins at around t = 100 kyr (models D4, D5, P098, P102, S090, S140, and S190). 
Recombination in the HII region is included in all cases. Here, we concentrate on the case 
with an imposed density perturbation (model D5) to illustrate why the perturbation grows 
when the amplitude of the imposed perturbation exceeds a critical value. 

Figures 5 (b)-(i) show a series of the 2D hydrogen number density plots (color) with 
incident photon flux contours (solid curves) at different times. The high density regions 
correspond to molecular gas compressed by the ablation pressure. A typical hydrogen number 
density in the HII region near the IF is 10^ cm"^ as shown in Figure 2. The solid white lines 
correspond to incident photon number contours from to 5 x 10^^ cm^^s^^ at intervals of 
1 X 10^^ cm^^s^^, starting from the ionization front where the ionization fraction (/) goes to 
zero. Since we do not show the whole computational domain, all the photon flux contours 
are not shown in the figures in the later phases. 



In Fig. 5 (b), the photon flux contours near to the IF are of the opposite phase as the IF 
because of the oscillation of the IF. This oscillation occurs by the stabilization mechanism 
seen in the linear regime (Paper I). In Fig. 5 (c), the | J| = photon flux contour separates 
slightly from the highly compressed molecular cloud around x = 0.23 pc, y = 2.3 pc (see 
also the close-up view in figure 6). This is the separation of the IF from the cloud surface. 
It happens because all the incident photons moving toward this region are absorbed by the 
recombined neutral hydrogen accumulated near the bubble vertex at a; = 0.23 pc. This could 
be understood as an extreme case of the normal stabilization mechanism of the IF instability, 
because the increasing density around the cavity causes strong absorption of the incident 
photons. During the separation of the IF, the cloud surface is shadowed. This shadowed 
region does not feel any ablation pressure because no absorption of the incident photons 
occur there, whereas the other parts of the cloud surface which are directly ablated by the 
incident photons feels strong ablation pressure. 

The region between the shadowed cloud surface and the separated IF becomes neutral. 
The appearance of the neutral region happens because the timescale for recombination is 
roughly (a^nj)"^ ~ 10^ yr for ionized hydrogen at a number density of n^ = lO^cm"^. Thus, 
once the cloud surface is shadowed, the neutral region quickly appears close to the shadowed 
surface. We can identify this region as a warm HI region (T ~ 4000 K ) because it has lost 
half of its thermal energy by recombination cooling. The pressure in this warm HI region 
is a few times smaller than that of the HII region and the molecular gas. This region has 
similar physical properties as the warm HI layer in front of PDRs described by Hollenbach 
& Tielens (1997) (see their Figure 3), which also arises due to a lowering of the ionization 
fraction. 

The appearance of this separation strongly affects the dynamics of the IF. In reaction, 
the compressed cloud expands into the low pressure region, namely, in the y direction, 
which triggers the second harmonic of the imposed initial perturbation in subsequent frames 
(Fig. 5(d) - (e)). Although the second harmonic of the imposed initial perturbation of the 
cloud surface can be clearly seen at t = 200 kyr, the second photon contour {\J\ = 1 x 
10^^ cm~^ s~^) from the IF still shows just a single mode. In other words, the density 
around the side (x = 0, 0.46 pc) is high, and the density around the center {x = 0.23) is low, 
although both regions on the cloud surface are spikes. Thus the stabilization mechanism 
seen in the linear regime (Paper I, and section 3.2) does not work in this case. The ablation 
pressure is strong around the center {x = 0.23 pc), because the two last contours are close to 
each other there (meaning the photon flux incident on the cloud surface is relatively high). 
As a result, a large amount of thermal energy is deposited there and compression by the 
strong ablation pressure occurs. By contrast, around the side {x = 0, 0.46 pc) the last two 
incident photon contours are separated from each other (meaning a lower photon flux reaches 
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the ablation front at the cloud surface). The resulting lower ablation pressure causes slight 
local expansion in the y direction (Fig.5 (e) and (f)). This expansion allows the growth of 
the spike at the sides (Fig.5 (g)), and a strong second harmonic has grown up. 

By 320 kyr, separation of the IF from the ablation front (cloud surface) near the lat- 
eral boundaries has occurred (Fig.5(h)), and by 400 kyr (Fig.5(i)), this separation has also 
appeared near the middle. The net result is that the perturbation amplitudes at the cen- 
ter and lateral edges grow. Once a critical amplitude is exceeded, at around 400 kyr, the 
stabilization mechanism is weakened by the steepness of the perturbations (i.e., the photon 
flux and therefore ablation rate along the steep sides of the perturbations drops, so there is 
much less density rise in the HII gas in the concave region), and the RT evolution proceeds 
rapidly into the deep nonlinear regime (Fig.5(j)). 

For the other cases with a large amplitude of initial perturbations (D4, P098, P102, S090, 
S140, and S190), we observe similar dynamics of the IF and cloud surface. The separation 
of the IF occurs, a second harmonic of the imposed perturbation is triggered, and finally two 
pillars grow. The center panels of Fig. 7 show a large perturbation growth in later phase of 
the dynamics using the different methods of the imposing the initial perturbation : density 
modulation (model D5; top), perturbation in photon number flux (model P102; middle), 
and surface perturbation (model S090; bottom). The growth of the second harmonic mode 
of the perturbation is observed in all cases. 



3.4. Velocity Gradient 

Figure 8 shows the t/-component of the velocity in the pillar for an initial perturbation 
in density (model D5), in photon flux contour (P102), and in the surface (S090). The Vy 
gradients are ~ 8.0 (D5), 16 (P102), and 12 (S090) km s~ pc~^. In a real pillar, an observer 
would see any velocity gradient corresponding to VySm.{i)^ where i is the inclination angle 
of the pillar with the plane of the sky. Indeed, Pound (1998) measured such gradients in 
the pillars of the Eagle Nebula. This observed velocities in Pillar II of the Eagle Nebula are 
also plotted for comparison. For the purposes of the flgure, we have taken i = 15 degrees, 
dividing the observed values by sin(15). Thus, for small inclination angles, our numerical 
results are good agreement with the observed Pillar II velocity gradient of 2.2 km s~^pc~^ 
(Pound (1998)). 
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3.5. Resolution Effects 

The calculation without the recombination effect, which was shown in Paper I, shows 
a large growth of the imposed perturbation. The surface also shows growth of a shorter 
wavelength perturbation due to purely numerical reasons, such as less dissipation. It is 
impossible to impose exactly one single mode using Cartesian coordinates; small amplitude 
and short wavelength perturbations will also always be imposed. The wavelength of the short- 
wavelength "noise" which grows depends on the resolution of the calculation. To eliminate 
the possibility of this effect in the current calculations, we also performed the calculations 
with higher and lower resolution in the x direction; specifically, twice (384) and a half (96) 
the number grid points than in our standard calculation. Figure 7 shows the results of higher 
(right panels) and lower (left panels) resolution calculations. As a reference, the results of 
the standard resolution calculation are also shown (center). The results of the three different 
ways for imposing perturbations are shown (top : t = 460 kyr, 50% peturbation in density; 
middle : t = 440 kyr, 30% peturbation in photon number flux during 102-110 kyr; bottom : 
t = 460 kyr, surface perturbation of amplitude 9 x 10~'^ pc). The pillars appear and grow at 
all resolutions, differing only in flne-scale structure at the cloud surfaces. The separation of 
the ionization front is also observed at the early phase of the dynamics in each case. Thus, 
we conclude that the pillar appearance is not a resolution-dependent effect. 



3.6. Discussion 

That we produce pillars is purely the result of our periodic boundary conditions. Unlike 
the pillars in our simulation, the Eagle pillars are not co-planar and do not share the same 
inclination angle (Pound (1998)). They are effectively three isolated pillars, rather than a 
single structure. Pillars in other HII regions are also typically isolated. However, our main 
aim in this paper has been to describe the instability of accelerating IF, so we have assumed 
initially a single wavelength and periodic boundary conditions. We have chosen only one 
set of free parameters — the wavelength and number density, initial magnetic pressure of the 
initial cloud, incident photon number flux at the boundary, and thickness of the cloud — 
although these speciflc values match those of Eagle Nebula reasonably well. 

The main key to this instability is the mismatch of perturbation modes between the IF 
and density in the HII region. The separation of the IF from the cloud is the trigger for 
this mismatch. Thus, it should still be possible to create an isolated pillar in the HII region 
under conditions where a single wavelength grows and becomes a pillar. Different wavelengths 
and multiwavelength cases should be investigated, however, since the phenomenon relies on 
nonlinear dynamics. 
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The effect of varying the magnetic field pressure sliould also be investigated, although 
the magnetic pressure is introduced just to prevent too high a number density from forming 
and does not affect the physical mechanism of the growth of the perturbation. An MHD 
calculation is one avenue for improvement over the current models. 

In this study we have adopted the on-the-spot approximation which ignores the recom- 
bination of the ionized hydrogen to the ground state, assuming that the photon emitted by 
this recombination is absorbed locally because its energy is high enough to ionize another hy- 
drogen atom. However, this approximation is not sufficient for Cartesian calculations when 
the directions of IF propagation and of the incident light rays become close to parallel. In 
such cases, irradiation by diffuse photons is the main source of ionization of the sides of the 
pillars. We expect to include this diffused photon transport in future calculations. 



4. CONCLUSION 

We have investigated hydrodynamic instability of the accelerating IF including detailed 
energy accounting. Three different ways of imposing the initial perturbation were adopted: 
(1) perturbation in density: (2) perturbation in photon number flux: and (3) surface per- 
turbation. 

We observe no large growth when the incident perturbation is small. This stabilization 
is the same as that observed in Paper I, and identical to the stabilization mechanism theoret- 
ically found by Axford (1964) for the non-accelerating IF. The strong absorption of incident 
photons causes decreasing ablation pressure around the cavity, which effectively smooths out 
small amplitude perturbations. 

In contrast, large growth of the second harmonic of the imposed perturbation is ob- 
served when the initial perturbation is large, namely, the ratio of the initial amplitude to 
wavelength is greater than 0.02. A slight separation of the IF from the cloud surface trig- 
gers the appearance of a warm HI region and the appearance of a second harmonic of the 
imposed perturbation. The separation of the IF from the cloud surface happens when all 
the incident photons traveling toward this region are absorbed by the recombined neutral 
hydrogen accumulated near the bubble vertex. The gas at this cloud surface at location 
begins to expand in the y direction during the separation because the pressure in the warm 
HI region is less than that in other regions. Although the second harmonic is imprinted 
on the cloud surface, the density perturbation and incident photon flux contours in the HII 
region are still single mode. The mismatch between modes at the cloud surface and in the 
density in the HII region above the cloud surface prevents stabilization and cause nonlinear 
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dynamics of the molecular cloud. In the later phase of the dynamics, large growth of the 
perturbation is observed. 

This new result for the IF instability does not depend on the method by which the 
initial perturbation is imposed. We also show that it is independent of the resolution of the 
calculation. The velocity gradient in the model pillar in the later phase of the dynamics is 
in good agreement with that observed in the Eagle Nebula. 

We would like to thank to Takashi Hosokawa, Mark Wolfire, Nobuhiko Izumi, and Robin 
Williams for useful discussions and comments. We acknowledge the anonymous referee for 
helpful comments which improved this manuscript. This work was performed under the 
auspices of the U.S. Department of Energy by the Lawrence Livermore National Laboratory 
under Contract No. W-7405-ENG-48 and with support from NASA Grant NRA 00-01-ATP- 
059 and from National Science Foundation under Grant No. AST-0228974. 
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Fig. 1. — Schematic diagram showing the position of the O stars that are photoevaporating 
the Eagle Nebula molecular cloud taken from Fig. 1 in Pound (1998). The size of the 
symbols indicates the relative strength of the Lyman continuum flux from each star. Most of 
the ionizing radiation comes from two stars of spectral type 05 V and 05.5 V in the nearby 
young cluster NGC 6611. 
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Fig. 2. — One dimensional hydrogen number density profiles of the case without perturba- 
tions at five different times, t = 0, 50, 100, 110, and 150 kyr. The profiles at three phases 
of the dynamics, compression (0 < t < 100 kyr), rarefaction (100 ^ ^ ^ HO kyr), and 
acceleration (t > 110 kyr) are shown. 
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Fig. 3. — One dimensional hydrogen number density profiles of the initial semi-infinite cloud. 
The cases with and without magnetic pressure at t = 500 kyr are shown. The number density 
of the compressed cloud without magnetic pressure is about 1.5 times higher than that with 
magnetic pressure. 
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Fig. 4. — The time evolution of the amphtude for different methods of imposing the initial 
perturbation. Top: the initial perturbation is imposed in density. Middle: the perturbation 
is imposed in the photon number flux. Bottom: a surface perturbation is imposed. Note 
the amplitude is defined as half the peak-to-valley height of the perturbation in the / = 
0.5 contour in the y direction. This does not always corresponds to the amplitude of the 
perturbation on the surface of the molecular cloud due to the "separation effect" (see text). 
When the imposed initial perturbation becomes large enough, the perturbations grow in all 
cases. When the initial perturbation is small, the growth is stabilized by recombination. 
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Fig. 5. — Log density (color) and incident photon flux contours (white sohd curves) from 
to I J| = 5 X 10^^ cm~^s~^ with uniform intervals of | J| = 1 x 10^^ cm~^s~^. The incident 
photon flux, which is parallel to the y axis, comes in from the top boundary. The dashed 
arrows in (b) indicate the initial photon flux from the top boundary, (a) A 30% initial 
perturbation in density is imposed (model D3), and no growth has occurred even at late 
times {t = 460 kyr). (b-j) Similar to (a) except than the initial density perturbation is 50% 
(model D5). The times corresponding to each 2D density snapshot are indicated at the top 
of each plot, (b) A shock is passing through the cloud, as indicated by the arrow. Note the 
\J\ = 1 X 10^^cm~^s~^ contour has the opposite phase as the IF. (c-g) The shock breaks out 
the back of the molecular cloud, and the acceleration phase begins. The solid arrows in (c) 
and (h) indicate the IF separation from the high density molecular cloud "surface", (h)-(i) 
Due to recombination, the IF separates from the high-density molecular cloud surface, which 
locally lowers the ablation pressure, and in (i) allows the perturbation spike in the center 
to grow due to the local pressure difference, (j) Once the perturbation growth exceeds a 
critical amplitude, the stabilization effect is reduced, and the evolution quickly enters the 
deep nonlinear regime, creating a pillar. 
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Fig. 6. — A close up of Fig. 5(c) around the IF. Density (x > 0) and pressure (x < 0) are 
plotted (color). The contours of the incident photon flux (|J| = 0, 1 x 10^^,2 x 10^^,3 x 
10^^,4 X 10^^ cm~^s~^) are plotted (x > 0). The contour of the ionization degree (/ = 0.5) 
is plotted (x < 0). The arrows in the right panel are velocity field. The IF slightly separated 
from the cloud surface at around x = 0.23, y = 2.3 pc, since all the incident photons are 
absorbed by the recombined hydrogen in the high density region (around x = 0.23, y = 2.4). 
A slight expansion of the cloud in the y direction can be seen. 
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Fig. 7. — Lower (left panels: half the number of grid points in x), standard zoning (middle 
panels) and higher (right panels: twice the number of grid points in x) resolution calculations. 
The contours show the same physical values as in Fig.5. The top panels are model D5 at 
t = 460 kyr, the middle panels are model 102 at t = 440 kyr, and the bottom panels are 
S090 at t = 460 kyr. All the cases show large growth of the second harmonic of the imposed 
initial perturbation, differing only in the fine-scale structure. 
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Fig. 8. — The velocity profile in the y direction at a; = 0.23 pc, for the different methods 
of initial perturbation: D5 (t = 460 kyr), P102 (t = 480 kyr), and S090 (t = 480 kyr). 
The pillar is oriented such that the right side represents its top and the left side represents 
its bottom. The points with error bars show the measured radial gas velocities in Pillar 
II of the Eagle Nebula from Pound (1998), divided by sin(15) to account for the line-of- 
sight inclination. The velocity gradients of models are in good agreement with the observed 
gradient. 



